1 Introduction

1.1 Objectifs du TP

En 2018, le PSAR analyse urbaine, ancêtre de la section analyse urbaine (DG/DDAR/DAR/DSAU), a développé un package R, nommé btb.

Sa principale fonction, kernelSmoothing, permet de réaliser très facilement un carroyage et un lissage sur des données géolocalisées. Grâce à différentes extensions, ce package peut être utilisé de plusieurs façons :

avec R
avec Python : libraryBTBpy
avec Qgis grâce à un plug-in développé par Lionel Cacheux (DR Grand Est)
avec ALICE, application intranet Insee développée par le PSAR Analyse urbaine (non maintenue)

À partir de données ponctuelles, nous allons apprendre :

À carroyer les informations.
À réaliser des lissages de densité, des lissages de moyennes, des lissages de taux et des lissages quantiles.
À calculer un indicateur sur une zone à façon à partir des données carroyées de l’Insee.

Les traitements sont réalisés avec le langage R. Le code de cette présentation est disponible ici.

1.2 Avertissements

1.2.1 Secret statistique

Avant toute diffusion auprès des partenaires, il faut bien veiller à respecter :

le secret
- primaire
- secondaire
- fiscal
les conventions établies avec les fournisseurs des données

1.2.2 Qualité des cartes

Pour simplifier les programmes présentés dans ce TP, les représentations graphiques ne respectent pas les règles élémentaires de la sémiologie cartographique très bien résumées par cette infographie :

Auteur : Timothée Giraud, auteur de la librairie mapsf

Il faudra bien veiller à appliquer ces règles de sémiologie sur les cartes que vous réaliserez ultérieurement.

1.2.3 Système de projection

Pour information, voici les systèmes de projection que vous pouvez régulièrement rencontrer pour la métropole :

Nom	Description	Code EPSG
Lambert93	Système de projection officiel pour la métropole	2154
LAEA	Système de projection européen	3035
WGS84	GPS (utile pour utiliser Leaflet)	4326

2 Configurations

2.1 Lancement de l’environnement de développement

Se connecter au SSP Cloud : trouver une doc toute faite quelque part
Lancer une session RStudio.
Faire un clone du projet git

2.2 Chargement des librairies

Quatre librairies sont nécessaires pour ce TP.

sf pour manipuler des fichiers spatiaux (importer des .shp, transformer des projections, et réaliser des géotraitements)
mapsf pour réaliser des cartes dans RStudio
leaflet pour réaliser des cartes interactives (fond de carte OpenStreetMap)
btb pour le carroyage et lissage
dplyr pour le traitement des données, en particulier l’agrégation géographique.

Remarque 2 : Le choix de dplyr plutôt que data.table se justifie ici du fait de sa forte compatibilité avec les objets géomatiques.

# Charger les librairies nécessaires

## Liste des librairies utilisées
packages <-  c("dplyr","sf","btb","mapsf","leaflet")

## Vérifier si la librairie est installée, si non l'installer, puis la charger
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

2.3 Chargement de la base

Pour ce TP, nous utilisons la base « Demandes de Valeurs Foncières », ou DVF, qui recense l’ensemble des ventes de biens fonciers réalisées au cours des cinq dernières années, en métropole et dans les départements et territoires d’outre-mer — sauf à Mayotte et en Alsace-Moselle. Les biens concernés peuvent être bâtis (appartement et maison) ou non bâtis (parcelles et exploitations). Les données sont produites par la direction générale des finances publiques. Elles proviennent des actes enregistrés chez les notaires et des informations contenues dans le cadastre.

A partir de cette source, nous avons constitué une base de données qui s’intéresse uniquement au périmètre de la petite couronne parisienne (départements 75, 92, 93 et 94) et au millésime 2021.

La base contient les 8 variables suivantes :

id_mutation : identifiant unique de la vente
date_mutation : date de la vente
type_local : appartement ou maison
nombre_pieces_principales : Nombre de pièces dans le logement
valeur_fonciere : prix de vente
surface_reelle_bati : surface du logement
x : longitude (en projection Lambert 93)
y : latitude (en projection Lambert 93)

Le code ayant permis de constituer la base de données est disponible ici : XXXXXX. TODO. + préciser éventuellement à grands traits les filtres réalisés dans la remarque ci-dessous.

Remarque importante : Cette base a été filtrée de manière à être la plus pédagogique possible pour cette formation. Elle n’est donc ni représentative de la réalité ni fidèle au champ proposé par le producteur. Mieux vaut donc de pas appuyer vos investissements immobiliers sur les résultats de ce TP.

Le code ci-dessous permet d’importer la première base de données ventesImmo_couronneParis.RDS utilisée dans ce tutoriel. Elle est stockée sous Minio, dans le “bucket public” : s3/kantunez/diffusion/projet_formation/r_lissage_spatial/.

# Charger la source de données (variable nommée dfBase) avec S3
url_bucket <- "https://minio.lab.sspcloud.fr/"
bucket <- "kantunez"
object <- "diffusion/projet_formation/r_lissage_spatial/ventesImmo_couronneParis.RDS"

dfBase <-
  aws.s3::s3read_using(
    FUN = base::readRDS,
    object = object,
    bucket = bucket
    ,
    opts = list("region" = "")
  )

Il est aussi possible de charger cette base en dehors du SSPCloud en la téléchargeant grâce à son URL puis en la chargeant dans R.

#Se positionner sur le répertoire de travail.
#Non nécessaire si projet RStudio créé au préalable
setwd("Mon/repertoire/qui/va/bien")

# Charger la source de données (variable nommée dfBase) avec son URL
download.file(paste0(url_bucket,bucket,"/",object), destfile = "ventesImmo_couronneParis.RDS")
dfBase <- readRDS("ventesImmo_couronneParis.RDS")

On peut ensuite manipuler notre base de données chargée.

# Visualiser les premières lignes de la base
head(dfBase)

##   id_mutation date_mutation  type_local nombre_pieces_principales valeur_fonciere
## 1 2021-447023    2021-01-08 Appartement                         3          480000
## 2 2021-447024    2021-01-05 Appartement                         2          345000
## 3 2021-447025    2021-01-08 Appartement                         3          384000
## 4 2021-447027    2021-01-08 Appartement                         3          261900
## 5 2021-447029    2021-01-05 Appartement                         2          407200
## 6 2021-447030    2021-01-15 Appartement                         3         1040000
##   surface_reelle_bati        x       y
## 1                  64 647357.3 6868635
## 2                  43 644483.5 6867695
## 3                  41 648001.8 6866153
## 4                  44 647414.8 6868937
## 5                  24 646929.9 6864730
## 6                  90 645132.7 6864646

dim(dfBase)

## [1] 34489     8

2.4 Chargement des données cartographiques

2.4.1 Fond de carte du territoire à étudier

Nous allons charger le contour géographique des départements de la petite couronne parisienne. En cartographie, ces fichiers sont appelés des « couches vectorielles » (à ne pas confondre avec les couches « raster » qui correspondent, pas exemple, aux fonds de carte satellite OpenStreetMap ou Google Maps). Une fois les départements chargés, nous allons sélectionner le contour de la commune de Paris.

Différents types de fichiers permettent de stocker des couches vectorielles. Les deux principaux sont le Geopackage (.gpkg) et le Shapefile (.shp). Nous recommandons l’utilisation du .gpkg pour les raisons suivantes :

Un .gpkg est un fichier unique alors qu’un fichier .shp tient dans 5 fichiers séparés et interdépendants.
Le format Géopackage est libre et ouvert, contrairement à Shapefile.
Un Shapefile impose des noms de variables de moins de 10 caractères, ce qui peut être pénible en pratique.
Le géopackage est devenu le format standard dans Qgis depuis la version 3 (logiciel libre de cartographie préconisé à l’Insee).

Pour aller plus loin sur la comparaison des formats, voir ici.

Récupérons le contour du territoire étudié avec la fonction read_sf du package sf.

# Charger le fond de carte du territoire étudié
object <- "diffusion/projet_formation/r_lissage_spatial/depCouronne.gpkg"
depCouronne_sf <- st_read(paste0(url_bucket,bucket,"/",object))

## Reading layer `depCouronne_sf' from data source 
##   `https://minio.lab.sspcloud.fr/kantunez/diffusion/projet_formation/r_lissage_spatial/depCouronne.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 4 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
## Projected CRS: RGF93 / Lambert-93

# Visualisons cette couche vectorielle
head(depCouronne_sf)

## Simple feature collection with 4 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
## Projected CRS: RGF93 / Lambert-93
##   code           libelle reg surf                           geom
## 1   75             Paris  11  105 MULTIPOLYGON (((660897 6860...
## 2   92    Hauts-de-Seine  11  176 MULTIPOLYGON (((648796 6847...
## 3   93 Seine-Saint-Denis  11  237 MULTIPOLYGON (((659428 6861...
## 4   94      Val-de-Marne  11  245 MULTIPOLYGON (((656908 6846...

plot(depCouronne_sf$geom)

depCouronne_sf <- depCouronne_sf %>% rename(geometry=geom)
# Sélection de Paris
paris_sf <- depCouronne_sf[depCouronne_sf$code=="75",]

# Visualisons cette couche vectorielle
head(paris_sf)

## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
## Projected CRS: RGF93 / Lambert-93
##   code libelle reg surf                       geometry
## 1   75   Paris  11  105 MULTIPOLYGON (((660897 6860...

plot(paris_sf$geometry)

Pour être certain que les données et le territoire soient dans le même système de projection, il est possible de transformer ce dernier à l’aide de la fonction st_transform.

# Transformer la projection
depCouronne_sf <- sf::st_transform(depCouronne_sf, crs = 2154)
paris_sf <- sf::st_transform(paris_sf, crs = 2154)

2.4.2 Territoires englobants, sélection des données à lisser

À partir de ces territoires, nous pouvons déterminer le rectangle englobant pour réduire la quantité de données manipulées.

Pour éviter les effets de bord, il faut également sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros). Autrement dit, si on ne sélectionne que les ventes immobilières situées dans Paris intramuros, et qu’on lisse ce nuage de points, les zones situées à proximité du périphérique seront artificiellement peu denses en ventes immobilières. En effet, elles sont situées à proximité de zones artificiellement vides en ventes immobilières dans notre nuage de points, à savoir dans les communes limitrophes. Pourtant, dans la réalité, il y a bien des ventes réalisées à Malakoff, Vanves, Montrouge, etc. Et ces ventes influencent la densité lissée aux niveau de la Porte de Vanves.

# Mise en forme de la couche "territoire" : sélection des variables et renommage
paris_sf$nom <- "territoire"
paris_sf <- paris_sf[,c("nom","geometry")]

# Création d'une bbox vectorielle autour du territoire
bbox <- sf::st_bbox(paris_sf)
bbox

##    xmin    ymin    xmax    ymax 
##  643076 6857499  660897 6867034

bbox_sf = st_sf(nom="bbox", geometry = st_as_sfc(bbox), crs = 2154)
bbox_sf

## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
## Projected CRS: RGF93 / Lambert-93
##    nom                       geometry
## 1 bbox POLYGON ((643076 6857499, 6...

# Création d'un buffer de la bbox, avec une marge de 2000 mètres

marge <- 2000
bufferBbox <- bbox
bufferBbox[["xmin"]] <- bufferBbox[["xmin"]]-marge
bufferBbox[["xmax"]] <- bufferBbox[["xmax"]]+marge
bufferBbox[["ymin"]] <- bufferBbox[["ymin"]]-marge
bufferBbox[["ymax"]] <- bufferBbox[["ymax"]]+marge
bufferBbox_sf = st_sf(nom="buffer_bbox", geometry = st_as_sfc(bufferBbox), crs = 2154)

# Assemblage des 3 couches dans une même table sf pour la cartographie
couches_plot <- rbind(bbox_sf,
      bufferBbox_sf,
      paris_sf)

# Passage d'une couche de polygones à une couche de lignes (pour la cartograhie)
couches_plot <- st_cast(couches_plot,"MULTILINESTRING")

# Cartographie pédagogique avec le package mapsf
mf_init(x = couches_plot, theme = "agolalight")
mf_map(x = st_cast(couches_plot[couches_plot$nom=="territoire",],"POLYGON"),
       type = "base", 
       col="wheat",
       border=FALSE,
       add = TRUE) 
mf_map(x = couches_plot,
       type = "typo", 
       var = "nom",
       pal = c("aquamarine4", "yellow3","wheat"),
       lwd=8,
       leg_title = "Zones d'analyse",
       add = TRUE) 
mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, IGN, mapsf")

Ce schéma présente :

en beige, le territoire d’intérêt (Paris)
en vert, la bbox du territoire d’intérêt (à savoir le plus petit rectangle englobant cette zone)
en jaune, la bbox élargie avec une marge de 2km. Autrement dit la zone tampon autour de la bbox permettant de prendre en compte des observations au-delà de la seule zone d’intérêt, et ainsi éviter les effets de bord au moment du lissage.

Pour sélectionner les données individuelles d’intérêt pour le lissage, on peut tout d’abord filtrer les observations dont les coordonnées sont comprises à l’intérieur du rectangle jaune. Ce filtre est très efficace computationnellement car il ne dépend que des valeurs numériques prises par les variables de longitude et de latitude et n’utilise pas les propriétés vectorielles des données qui sont des attributs plus chronophages à utiliser.

# Repérer les coordonnées extrêmes, avec une marge (ici 2km)
bufferBbox

##    xmin    ymin    xmax    ymax 
##  641076 6855499  662897 6869034

xMin <- bufferBbox["xmin"]
xMax <- bufferBbox["xmax"]
yMin <- bufferBbox["ymin"]
yMax <- bufferBbox["ymax"]

# Ne garder que les données dans le rectangle englobant le territoire étudié, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= xMin & dfBase$x <= xMax & dfBase$y >= yMin & dfBase$y <= yMax, ]
dim(dfBase_filtre)

## [1] 24419     8

Il est possible de procéder autrement. On peut utiliser nos données individuelles comme un ensemble de points (approche vectorielle) et procéder à des intersections géographiques. Pour cela :

on créer une zone tampon (buffer) autour du territoire d’intérêt, avec une marge (ici 2 000m), sous la forme d’un objet sf vectoriel
On repère les observations comprises dans cette zone tampon par intersection géographique, après avoir transformé nos observation en points vectoriels.

Cette méthode, bien qu’élégante, peut-être lourdre d’un point de vu computationnel. Avec des données volumineuses, il convient au minimum de faire un premier filtrage avec la précédente méthode.

# Transformer le tableau de données filtré en objet géographique
sfBase <- sf::st_as_sf(dfBase, coords = c("x", "y"), crs = 2154)

# Création d'un buffer autour du territoire
buffer_sf <- st_buffer(paris_sf, dist = 2000)

# Repérer les indices des observations contenues dans notre buffer d'intérêt
indiceObsContenues <- unlist(sf::st_contains(buffer_sf, sfBase))

# Réduire la base aux seules observations dans le territoire
dfBase_filtre <- dfBase[indiceObsContenues, ]

# Mise en forme de la couche buffer
buffer_sf$nom <- "buffer"

# Assemblage des 2 couches dans une même table sf pour la cartographie
couches_plot <- rbind(buffer_sf,paris_sf)

# Passage d'une couche de polygones à une couche de lignes (pour la cartographie)
couches_plot <- st_cast(couches_plot, "MULTILINESTRING")

# Cartographie pédagogique avec le package mapsf
mf_init(x = couches_plot, theme = "agolalight")
mf_map(x = couches_plot,
       type = "typo", 
       var = "nom",
       pal=c("yellow3","wheat"),
       lwd=8,
       leg_title = "Zones d'analyse",
       add = TRUE) 
mf_map(x = st_cast(couches_plot[couches_plot$nom=="territoire",],"POLYGON"),
       type = "base", 
       col="wheat",
       border=FALSE,
       add = TRUE) 

# Échantillon de 500 observations dans le buffer
sfBase_sample <- sfBase[indiceObsContenues, ]
sfBase_sample <- sfBase_sample[sample(1:nrow(sfBase_sample),2000) ,]
mf_map(sfBase_sample,
       add = TRUE,
       leg_title = "observations")

mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque: Pour la zone tampon, il convient de prendre une marge légèrement plus grande que le rayon de lissage envisagé, ceci afin d’éviter les effets de bord tout en limitant les temps de calcul. Une fois le lissage réalisé, on peut réduire la carte à la zone étudiée.

3 Carroyage de données

Avant de lisser les données ponctuelles, on peut souhaiter représenter ces données sous forme carroyée afin de se les approprier. Il faut pour cela :

Tracer une grille carroyée du territoire en s’appuyant sur la répartition spatiale des points ;
Agréger les données à l’échelle des carreaux créés précédemment ;
Cartographier ces carreaux sous forme d’une carte choroplète (discrétisation des valeurs prises par chaque polygone, à savoir chaque carreau dans le cas présent).

Le code ci-dessous procède de la sorte :

on associe chaque point (= vente géolocalisée) au centroïde du carreau auquel il appartient (le territoire est découpé en carreaux de 200 mètres à partir de l’origine du référentiel)
On agrège les données sur chaque centroïde de la grille
On passe d’une table de centroïdes à une table de carreaux vectoriels grâce à la fonction dfToGrid du package btb.

La fonction dfToGrid attend comme paramètres obligatoires :

df : un tableau avec les colonnes x et y représentant les coordonnées des centroïdes de la grille. (Remarque : la grille n’est pas nécessairement régulière.)
sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé.

# Arrondir les coordonnées des observations aux centroides les plus proches
iCellSize = 200 # Carreaux de 200m
points_carroyage <- dfBase_filtre
points_carroyage$x_centroide = as.integer(floor(points_carroyage$x / iCellSize) * iCellSize + (iCellSize / 2))
points_carroyage$y_centroide = as.integer(floor(points_carroyage$y / iCellSize) * iCellSize + (iCellSize / 2))

head(points_carroyage)

##    id_mutation date_mutation  type_local nombre_pieces_principales valeur_fonciere
## 3  2021-447025    2021-01-08 Appartement                         3          384000
## 5  2021-447029    2021-01-05 Appartement                         2          407200
## 6  2021-447030    2021-01-15 Appartement                         3         1040000
## 9  2021-447034    2021-01-18 Appartement                         4          625000
## 10 2021-447035    2021-01-18 Appartement                         4          721250
## 11 2021-447037    2021-01-07 Appartement                         4         1035000
##    surface_reelle_bati        x       y x_centroide y_centroide
## 3                   41 648001.8 6866153      648100     6866100
## 5                   24 646929.9 6864730      646900     6864700
## 6                   90 645132.7 6864646      645100     6864700
## 9                   68 647494.6 6866019      647500     6866100
## 10                  94 643891.9 6864587      643900     6864500
## 11                  94 648239.0 6866235      648300     6866300

# Compter le nombre de ventes par carreau (si on cherche à représenter la densité)
points_carroyage <- points_carroyage %>% 
  group_by(x=x_centroide,y=y_centroide) %>% 
  count(name = "nbVentes")
  

# Générer la grille
points_carroyage <- btb::dfToGrid(df = points_carroyage, sEPSG = "2154", iCellSize = iCellSize)

# Restriction du champ : on ne retient que les carreaux intersectant Paris
points_carroyage <- points_carroyage[unlist(st_intersects(paris_sf,points_carroyage)),]

# Cartographie

mf_init(x=points_carroyage,theme = "agolalight")
mf_map(x = points_carroyage,
       type = "choro",
       var="nbVentes",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       leg_val_rnd = 1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Carroyage du nombre de ventes",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque 1 : Le carroyage peut aussi être utilisé pour simplifier les données avant lissage si le calcul de lissage est trop long en raison d’un très grand nombre d’observations.

Remarque 2 : Les étapes avant la génération de la grille ne sont pas intégrées dans des fonctions du package btb, il faut donc s’inspirer du code ci-dessus pour reproduire le même type de carroyage. Il est envisagé qu’une future version du package btb simplifie le code permettant de réaliser un carroyage.

4 Lissage

4.1 Calcul de densité

Un lissage simple à réaliser correspond au calcul d’une densité, à savoir une quantité par unité de surface (un carreau). Dans cette section nous allons calculer la densité ventes de logements sur la ville de Paris au cours de l’année 2021.

Pour ce faire, nous allons au préalable créer une variable nbObsLisse qui vaudra 1 pour chaque observation (donc logement vendu en 2021) et qui sera lissée avec la fonction kernelSmoothing du package btb.

Toujours pour éviter les effets de bord, le lissage sera effectué sur une zone plus large que la zone d’intérêt (utilisation d’un buffer).

4.1.1 Grille automatique; carreaux 200m; rayon de lissage 400m

Pour ce premier lissage, nous choisissons des carreaux de 200 mètres et un rayon de lissage de 400 mètres.

La fonction de lissage btb::kernelSmoothing permet de générer automatiquement le carroyage préalable au lissage.

# Créer une nouvelle variable pour compter le nombre de ménages lissé
dfBase_filtre$nbObsLisse <- 1

# Visualiser les premières lignes de la base
head(dfBase)

##   id_mutation date_mutation  type_local nombre_pieces_principales valeur_fonciere
## 1 2021-447023    2021-01-08 Appartement                         3          480000
## 2 2021-447024    2021-01-05 Appartement                         2          345000
## 3 2021-447025    2021-01-08 Appartement                         3          384000
## 4 2021-447027    2021-01-08 Appartement                         3          261900
## 5 2021-447029    2021-01-05 Appartement                         2          407200
## 6 2021-447030    2021-01-15 Appartement                         3         1040000
##   surface_reelle_bati        x       y
## 1                  64 647357.3 6868635
## 2                  43 644483.5 6867695
## 3                  41 648001.8 6866153
## 4                  44 647414.8 6868937
## 5                  24 646929.9 6864730
## 6                  90 645132.7 6864646

Pour réaliser le lissage, il suffit de fournir à la fonction kernelSmoothing quatre paramètres :

dfObservations : le tableau des données à lisser. Il doit nécessairement contenir une colonne x, une colonne y, et 1 à n colonnes numériques (variables à lisser) ;
sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé ;
iCellSize : un entier indiquant la taille des carreaux ;
iBandwidth : un entier indiquant le rayon de lissage.

Attention, la fonction retourne une erreur en cas de :

présence d’une variable non-numérique ;
valeur(s) absente(s) dans les colonnes x ou y.

# Lisser
dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 400)

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 650200 ymin: 6855400 xmax: 653400 ymax: 6855800
## Projected CRS: RGF93 / Lambert-93
##        x       y nbObsLisse                       geometry
## 1 651500 6855500  1.0477469 POLYGON ((651400 6855600, 6...
## 2 652100 6855500  0.8521979 POLYGON ((652000 6855600, 6...
## 3 652700 6855500  2.7218839 POLYGON ((652600 6855600, 6...
## 4 653300 6855500  1.4556305 POLYGON ((653200 6855600, 6...
## 5 650300 6855700  2.1375004 POLYGON ((650200 6855800, 6...
## 6 650700 6855700  2.0552302 POLYGON ((650600 6855800, 6...

# Nombre de lignes lissées
nrow(dfBaseLisse)

## [1] 3259

La fonction kernelSmoothing génère une grille de carreaux, et chaque carreau est affecté de la densité lissée en son centroïde. À noter que la grille carroyée épouse globalement le périmètre des points en entrée de la fonction (elle va même un peu au-delà par défaut, cf. séquence théorique sur le lissage de la formation “Comment utiliser les outils de l’analyse urbaine ?”).

mf_init(x=paris_sf,theme = "agolalight")
mf_map(x = st_cast(dfBaseLisse[,c("geometry")],"LINESTRING"), 
       lwd=1,add=T)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=8,
       col="wheat",add = TRUE)
mf_layout(title = "Gestion du territoire d'analyse",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Une propriété particulièrement importante de la fonction de lissage classique du package btb est qu’elle est conservative. Cela signifie que nous avons la même somme des variables additives sur le champ géographique concerné avant et après lissage.

sum(dfBase_filtre$nbObsLisse) # Nombre de ventes dans la petite couronne

## [1] 22221

sum(dfBaseLisse$nbObsLisse) # Après lissage, nombre lissé de ventes dans les carreaux

## [1] 22221

Remarque : Veillez toujours à respecter le secret statistique au moment de la publication de vos cartes ! Par exemple, si vous utilisez des données issues des bases fiscales, vous ne pouvez représenter des informations sur des carreaux comportant moins de 11 observations, même si ce nombre est lissé… Pour ce faire, vous pouvez réaliser le filtre ci-dessous :

# Exclure les carreaux ne respectant pas le secret statistique
dfBaseLisse <- dfBaseLisse[dfBaseLisse$nbObsLisse >= 11, ]

Voyons maintenant le nombre lissé de ventes de logements à Paris en 2021 suite à ce premier lissage.

Quelques informations préalables :

On cherche à analyser le phénomène sur le périmètre de la ville de Paris ;
Ainsi, on lisse les données de ventes sur un périmètre plus large que la ville de Paris afin d’éviter les effets de bord aux frontière de la commune. En effet, la base dfBase_filtre contient toutes les ventes comprises dans un buffer de 2 km autour de Paris.
Suite au lissage, une grille carroyée est produite comportant les données lissées. Les carreaux vont bien au-delà de la commune de Paris. Néanmoins, pour la représentation, on ne sélectionne que les carreaux intersectant Paris (grâce à la fonction sf::st_intersects).
Pour cartographier les carreaux :
- On réalise une carte de type choroplèthe ;
- Avec une méthode de discrétisation de la densité lissée : ici, par quantiles par exemple ;
- Avec 5 classes (c’est-à-dire 5 couleurs).

# Filtrage des carreaux lissés intersectant la ville de Paris
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]

# Carte lissée
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 400m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque : la variable lissée nbObsLiss correspond au nombre de ventes (observations) par carreau (ici de 200 mètres). Ainsi, il est possible d’obtenir une densité au km² en multipliant la variable nbObsLiss par 25 dans le cas présent.

4.1.2 Grille automatique; carreaux 200; rayon de lissage 600m

Il n’est pas possible de connaitre a priori la taille optimale du rayon de lissage. Il faut donc en essayer plusieurs pour trouver le meilleur compromis entre précision et généralisation. Nous allons donc faire le même lissage que précédemment avec un rayon de 600 m au lieu de 400 m.

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 600)

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
##        x       y nbObsLisse                       geometry
## 1 651300 6855300  0.1065710 POLYGON ((651200 6855400, 6...
## 2 651500 6855300  0.1384200 POLYGON ((651400 6855400, 6...
## 3 651700 6855300  0.2001949 POLYGON ((651600 6855400, 6...
## 4 651900 6855300  0.2094593 POLYGON ((651800 6855400, 6...
## 5 652100 6855300  0.1441657 POLYGON ((652000 6855400, 6...
## 6 652300 6855300  0.1852248 POLYGON ((652200 6855400, 6...

# Filtrage des carreaux lissés dans Paris
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]

# Carte lissée
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 600m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

À noter qu’on peut réitérer l’opération avec un rayon de lissage de 1000 mètres.

## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651000 ymin: 6855000 xmax: 652200 ymax: 6855200
## Projected CRS: RGF93 / Lambert-93
##        x       y nbObsLisse                       geometry
## 1 651100 6855100 0.09890581 POLYGON ((651000 6855200, 6...
## 2 651300 6855100 0.09483879 POLYGON ((651200 6855200, 6...
## 3 651500 6855100 0.11358059 POLYGON ((651400 6855200, 6...
## 4 651700 6855100 0.14485863 POLYGON ((651600 6855200, 6...
## 5 651900 6855100 0.16973038 POLYGON ((651800 6855200, 6...
## 6 652100 6855100 0.18790912 POLYGON ((652000 6855200, 6...

À mesure que nous augmentons le rayon de lissage, la carte révèle des aspects structurels des données. Mais ceci se fait au détriment des spécificités locales visibles seulement avec un petit rayon de lissage. Le choix revient dont au statisticien-géographe, en fonction de sa connaissance des données et des objectifs recherchés.

Enfin, la grille carroyée est ici visible à des fins pédagogiques, il est bien-entendu possible et conseillé de la supprimer sur les cartes publiées comme dans cet exemple :

mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       border = NA, # C'est ici que ça se passe
       # lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m, sans grille",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

4.1.3 Exemple d’effet de bord à éviter

Ci-dessous, nous illustrons les effets de bord produits si on lisse uniquement les ventes réalisées dans Paris intramuros. Dans les deux cartes ci-dessous, le rayon de lissage (2 000 mètres) et les bornes de discrétisation sont communs pour assurer la comparabilité des deux cartes.

Dans le premier cas, on lisse les ventes situées dans Paris et ses alentours. - Dans le second cas, on ne lisse que les ventes réalisées dans intramuros.

Ainsi, on peut remarquer que les effets de bord se matérialisent à la périphérie de Paris, où les densités lissées sont artificiellement plus faibles que quand on lisse en prenant une marge.

# Lissage avec gestion des effets de bord
dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 2000)


# Lissage sans gestion des effets de bord (Base des ventes uniquement dans intramuros)
sfBase_intramuros <- sfBase[unlist(st_contains(paris_sf,sfBase)),]
dfBase_intramuros <- sfBase_intramuros %>% 
  mutate(x=st_coordinates(geometry)[,1],
         y=st_coordinates(geometry)[,2],
         nbObsLisse=1) %>% 
  st_drop_geometry() %>% 
  select(nbObsLisse,x,y)

dfBaseLisse_intramuros <- btb::kernelSmoothing(dfObservations = dfBase_intramuros, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 2000)

Évoquer le paramétrage de la grille carroyée afin de prendre en compte les limites artificielles (limites administratives) et les limites physiques (mer, montagne, etc.)

4.2 Calcul de moyenne

Dans cette section nous allons nous intéresser au prix moyen des logements vendus à Paris en 2021.

Remarque : Une moyenne n’est pas le lissage du rapport (liss(variable / nbObsLisse)) mais le rapport des lissages : liss(variable) / liss(nbObs).

Ainsi, il convient de :

Lisser les prix de ventes des biens vendus
Lisser le nombre de biens vendus (exactement comme précédemment)
Faire le ratio des deux précédents lissages pour chaque carreau de la grille produite.

Dans l’exemple ci-dessous, on prend un rayon de lissage de 600 m :

dataLissage <- dfBase_filtre[,c("nbObsLisse","valeur_fonciere","x","y")] # Ajout de la variable "valeur_fonciere"
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 600)

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
##        x       y nbObsLisse valeur_fonciere                       geometry
## 1 651300 6855300  0.1065710        38517.44 POLYGON ((651200 6855400, 6...
## 2 651500 6855300  0.1384200        47025.68 POLYGON ((651400 6855400, 6...
## 3 651700 6855300  0.2001949        66234.22 POLYGON ((651600 6855400, 6...
## 4 651900 6855300  0.2094593        72332.24 POLYGON ((651800 6855400, 6...
## 5 652100 6855300  0.1441657        52248.55 POLYGON ((652000 6855400, 6...
## 6 652300 6855300  0.1852248        68258.16 POLYGON ((652200 6855400, 6...

# Calculer le taux
dfBaseLisse$prixMoyen <- dfBaseLisse$valeur_fonciere / dfBaseLisse$nbObsLisse

# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris, 
       type = "choro",
       var="prixMoyen",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage des prix de vente avec un rayon de 600m",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

On remarque des prix moyens particulièrement élevés dans le centre et l’Ouest de la capital, ce qui semble conforme à l’intuition.

Attention cependant aux éventuelles valeurs atypiques (ventes extrêmement chères, erreur dans les données) : le lissage en moyenne est potentiellement sensible à ces valeurs atypiques (contrairement au lissage quantile, voir infra).

4.3 Calcul de taux

Dans cette section, nous allons nous intéresser à la proportion de ventes portant sur des grands logements (plus de 4 pièces principales).

# Créer une nouvelle variable pour détecter les logements de plus de 4 pièces principales
dfBase_filtre$quatrePieces <- ifelse(dfBase_filtre$nombre_pieces_principales >= 4, 1, 0)

dataLissage <- dfBase_filtre[,c("nbObsLisse","quatrePieces","x","y")]
dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 600)

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
##        x       y nbObsLisse quatrePieces                       geometry
## 1 651300 6855300  0.1065710   0.07578049 POLYGON ((651200 6855400, 6...
## 2 651500 6855300  0.1384200   0.07838126 POLYGON ((651400 6855400, 6...
## 3 651700 6855300  0.2001949   0.05794216 POLYGON ((651600 6855400, 6...
## 4 651900 6855300  0.2094593   0.02832061 POLYGON ((651800 6855400, 6...
## 5 652100 6855300  0.1441657   0.02332406 POLYGON ((652000 6855400, 6...
## 6 652300 6855300  0.1852248   0.13111821 POLYGON ((652200 6855400, 6...

# Calculer le taux
dfBaseLisse$txQuatrePieces <- dfBaseLisse$quatrePieces / dfBaseLisse$nbObsLisse

dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]

# Cartographie
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris, 
       type = "choro",
       var="txQuatrePieces",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"), 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Logements avec 4 pièces ou plus (rayon de 600m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

La carte des prix moyens et la carte des grands logements se ressemblent, fort logiquement. On peut alors avoir envie de lisser le prix au mètre-carré.

Pour ce faire, n’oubliez pas qu’il convient de lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur) des logements vendus. Autrement, si on lisse le prix au m² de chaque logement vendu, on sur-pondère artificiellement les prix au m² des petits logements (car l’unité statistique devient alors le logement, et non le m² vendu).

dataLissage <- dfBase_filtre[,c("surface_reelle_bati","valeur_fonciere","x","y")]

dfBaseLisse <- btb::kernelSmoothing(dfObservations = dataLissage,
                                    sEPSG = "2154",
                                    iCellSize = 200,
                                    iBandwidth = 600)

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 651200 ymin: 6855200 xmax: 652400 ymax: 6855400
## Projected CRS: RGF93 / Lambert-93
##        x       y surface_reelle_bati valeur_fonciere                       geometry
## 1 651300 6855300            8.067281        38517.44 POLYGON ((651200 6855400, 6...
## 2 651500 6855300           10.163131        47025.68 POLYGON ((651400 6855400, 6...
## 3 651700 6855300           13.385886        66234.22 POLYGON ((651600 6855400, 6...
## 4 651900 6855300           12.786188        72332.24 POLYGON ((651800 6855400, 6...
## 5 652100 6855300            7.943338        52248.55 POLYGON ((652000 6855400, 6...
## 6 652300 6855300           12.310115        68258.16 POLYGON ((652200 6855400, 6...

# Calculer le prix au m²
dfBaseLisse$prixM2 <- dfBaseLisse$valeur_fonciere / dfBaseLisse$surface_reelle_bati

# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
       type = "choro",
       var="prixM2",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Prix au m² (rayon de 600m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

4.4 Calcul de quantiles géographiquement pondérés

Dans cette section, nous allons calculer le 1er décile, la médiane et le 9ème décile du prix de vente des logements à Paris. Ce calcul se fait toujours avec la fonction kernelSmoothing mais il faut ajouter un paramètre, vQuantiles, qui contient la liste des quantiles souhaités.

Remarque : Le calcul de quantiles géographiquement pondérés présente certaines différences par rapport au lissage classique vu jusqu’à présent :

il n’est pas conservatif
la fonction crée
- une colonne nbObs indiquant le nombre d’observations réel (non lissé) ayant contribué au calcul du carreau.
- pour chaque variable, autant de colonnes qu’il y a de quantiles souhaités

On représente ici le lissage du premier décile et de la médiane.

Pour rappel, le lissage quantile (médiane, déciles . . . ) permet d’avoir des indicateurs moins sensibles aux valeurs extrêmes et d’enrichir l’analyse avec des indicateurs de dispersion (comme par exemple l’écart interquantile).

# Visualiser les premières lignes des données lissées
head(dfBaseLisse)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 650800 ymin: 6854800 xmax: 652000 ymax: 6855000
## Projected CRS: RGF93 / Lambert-93
##   nbObs valeur_fonciere_01 valeur_fonciere_05 valeur_fonciere_09      x       y
## 1    55             152000             335000             700000 650900 6854900
## 2    56             152000             319000             675000 651100 6854900
## 3    58             152000             319000             655000 651300 6854900
## 4    60             152000             319000             655000 651500 6854900
## 5    58             182000             324000             655000 651700 6854900
## 6    60             185000             324000             655000 651900 6854900
##                         geometry
## 1 POLYGON ((650800 6855000, 6...
## 2 POLYGON ((651000 6855000, 6...
## 3 POLYGON ((651200 6855000, 6...
## 4 POLYGON ((651400 6855000, 6...
## 5 POLYGON ((651600 6855000, 6...
## 6 POLYGON ((651800 6855000, 6...

# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
       type = "choro",
       var="valeur_fonciere_01",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Premier décile des prix de vente (rayon de 1500m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

# Cartographie
dfBaseLisse_paris <- dfBaseLisse[unlist(st_intersects(paris_sf,dfBaseLisse)),]
mf_init(x=dfBaseLisse_paris,theme = "agolalight")
mf_map(x = dfBaseLisse_paris,
       type = "choro",
       var="valeur_fonciere_05",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Médiane des prix de vente (rayon de 1500m)",credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Pour que le quantile lissé ait du sens, il faut qu’un nombre conséquent d’observations ait participé à sa construction. Ainsi, il est conseillé d’utiliser un rayon de lissage suffisamment élevé.

summary(dfBaseLisse_paris$nbObs)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      21     715    1029    1096    1430    2907

# Part de centroïdes dans le quantile lissé a été calculé avec moins de 10 observations
dfBaseLisse_paris %>% st_drop_geometry() %>% group_by(nbObs<10) %>% count()

## # A tibble: 1 × 2
## # Groups:   nbObs < 10 [1]
##   `nbObs < 10`     n
##   <lgl>        <int>
## 1 FALSE         2672

5 Indicateurs sur zonage à façon

Dans cette section nous allons produire des indicateurs sur des zones à façon.

Pour commencer, nous calculerons le nombre et le prix moyen des transactions immobilières en 2021 dans le “triangle d’or” de la capitale. Ensuite, nous mobiliserons les données carroyées à 200 mètres disponibles sur insee.fr) pour calculer la proportion de logements sociaux autour du White.

Remarque : une fonction généralisant ce second exercice est disponible en ligne, sur insee.fr, en accompagnement des données carroyées.

La zone à façon que nous considérons est le cercle de 1 km de rayon autour du White, à Montrouge. L’objectif est de déterminer la part de logements sociaux dans cette zone, à partir des données carroyées (carreaux de 200 mètres).

5.1 Indicateurs sur zones à façon à partir de données ponctuelles

Dans un premier temps, on importe la zone à façon. Ici, il s’agit du quartier appelé «~Triangle d’or~», situé dans le 8e arrondissement de Paris. Ce quartier se caractérise par ses immeubles cossus, ses hôtels et ses commerces de luxe.

object <- "diffusion/projet_formation/r_lissage_spatial/triangle_or.kml"
triangleOr <- st_read(paste0(url_bucket,bucket,"/",object))

## Reading layer `triangle_or' from data source 
##   `https://minio.lab.sspcloud.fr/kantunez/diffusion/projet_formation/r_lissage_spatial/triangle_or.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 2.300717 ymin: 48.86447 xmax: 2.310474 ymax: 48.87203
## Geodetic CRS:  WGS 84

leaflet(data=triangleOr) %>% addTiles() %>% addPolygons()

Remarque : ce fichier géogrpahique a pour extension .kml. Ce format se lit sans problème avec la fonction sf::st_read. Pour information, ce triangle a été crée manuellement sur le site du Géoportail.

Pour connaître le nombre et le prix moyen des transactions immobilières réalisées en2021 dans ce quartier, nous procédons par intersection géographique entre des points (les transactions) et un polygone (le quartier). Attention à bien s’assurer que les deux couches géographiques soient bien projetées dans le même système de projection.

# Système de projection de triangleOr
st_crs(triangleOr)[["input"]]

## [1] "WGS 84"

# Transformation du système de projection de triangleOr 
triangleOr <- st_transform(triangleOr,crs = 2154)

# Intersection géographique avec les transactions
pos_intersec <- st_intersects(triangleOr,sfBase) %>% unlist()
transac_triangle <- sfBase[pos_intersec,]

# Visualisation de points retenus
transac_triangle

## Simple feature collection with 20 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 648767.5 ymin: 6863055 xmax: 649194 ymax: 6863601
## Projected CRS: RGF93 / Lambert-93
## First 10 features:
##       id_mutation date_mutation  type_local nombre_pieces_principales valeur_fonciere
## 18629 2021-489256    2021-01-15 Appartement                         5         3300000
## 18676 2021-489336    2021-01-25 Appartement                         3          590000
## 18677 2021-489337    2021-01-28 Appartement                         4         1511345
## 18698 2021-489369    2021-02-05 Appartement                         6         2450000
## 18723 2021-489406    2021-01-20 Appartement                         3          997000
## 18725 2021-489410    2021-02-12 Appartement                         3         1605000
## 18748 2021-489452    2021-02-24 Appartement                         2          821500
## 18757 2021-489464    2021-02-04 Appartement                         2         1200000
## 18768 2021-489483    2021-02-19 Appartement                         2           74900
## 18777 2021-489495    2021-02-19 Appartement                         4         2500000
##       surface_reelle_bati                 geometry
## 18629                 180 POINT (648918.9 6863145)
## 18676                  53 POINT (648903.5 6863601)
## 18677                  65   POINT (649194 6863377)
## 18698                 133 POINT (648814.5 6863055)
## 18723                  72 POINT (648870.5 6863135)
## 18725                 100 POINT (648767.5 6863434)
## 18748                  60 POINT (648771.3 6863503)
## 18757                  74 POINT (648767.5 6863434)
## 18768                  44 POINT (648814.2 6863144)
## 18777                 100 POINT (648970.8 6863181)

# Cartographie des points dans le triangle
leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = triangleOr %>% st_transform(4326)) %>% 
  addMarkers(data = transac_triangle %>% st_transform(4326))

Remarque 1 : pour tracer cette dernière carte, nous devons re-projeter le triangle et les points en WGS84 (4326), car la fonction Leaflet ne le fait pas automatiquement (contrairement à mapview).

Remarque 2 : on ne voit que 12 points sur la cartes, alors qu’on a trouvé 20 transactions dans le triangle : ceci s’explique car certaines transactions sont géolocalisées au même endroit (appartements situés dans un même immeuble par exemple).

L’ensemble des transactions ayant été récupérées par intersections géographiques, nous pouvons enfin calculer les indicateurs souhaités.

# Nombre de transactions en 2021 : 
nrow(transac_triangle)

## [1] 20

# Prix moyen
mean(transac_triangle$valeur_fonciere)

## [1] 1637256

5.2 Indicateurs sur zones à façon à partir de données carroyées

Le principe de la démarche consiste à : 1. Déterminer, pour une zone donnée, l’ensemble des carreaux qui la recouvrent 2. Calculer les agrégats sur cet ensemble de carreaux.

Remarque : L’ensemble de carreaux étant en général plus large que la zone, les agrégats obtenus seront des estimations qui tendent à surestimer les valeurs réelles. Le problème ne se pose pas si vous travaillez directement sur des données disponibles au niveau “individu statistique” (“avec des x,y”).

@Julien : on peut aussi donner un exemple d’agrégation aux x,y avec les données de DVF juste avant l’exemple au carreau. Ca me paraitrait intéressant.

# library(stringr, warn.conflicts = FALSE)

# Charger les données carroyées de Filosofi 2015 à 200m
# Uniquement en région parisienne
# En ne conservant que 6 variables dont le nombre de logements sociaux
st_read_maison <- function(chemin_tab){
  requete <- "SELECT IdINSPIRE,Depcom,I_est_cr,Men, Log_soc, geom
            FROM Filosofi2015_carreaux_200m_metropole
            WHERE SUBSTR(Depcom, 1, 2) IN ('75','92','93','94') "
  sf::st_read(chemin_tab, query = requete)
}

object = "diffusion/projet_formation/r_lissage_spatial/Filosofi2015_carreaux_200m_metropole.gpkg"
carreaux <-  st_read_maison(paste0(url_bucket,bucket,"/",object))

head(carreaux)

## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 652623.2 ymin: 6858353 xmax: 655268.9 ymax: 6866380
## Projected CRS: RGF93 / Lambert-93
##                        IdINSPIRE Depcom I_est_cr Men Log_soc
## 1 CRS3035RES200mN2893400E3763200  75119        0 990     937
## 2 CRS3035RES200mN2890000E3762200  75111        0 926      80
## 3 CRS3035RES200mN2893400E3762400  75119        0 508     323
## 4 CRS3035RES200mN2891800E3763600  75119        0 633       0
## 5 CRS3035RES200mN2890800E3763200  75120        0 349     188
## 6 CRS3035RES200mN2885800E3760600  75113        0 751     344
##                             geom
## 1 MULTIPOLYGON (((654521.9 68...
## 2 MULTIPOLYGON (((653843.4 68...
## 3 MULTIPOLYGON (((653725.1 68...
## 4 MULTIPOLYGON (((655069.7 68...
## 5 MULTIPOLYGON (((654764.7 68...
## 6 MULTIPOLYGON (((652641.8 68...

# Création du cercle de 1km autour du White
white <- st_sfc(st_point(c(649218.36,6857569.14)),crs=2154)
buffer_white <- st_buffer(white,1000)

# Sélection des carreaux intersectant notre zone d'intérêt
carreaux_select <- carreaux[unlist(st_intersects(buffer_white,carreaux)),]

# Cartographie

leaflet() %>% 
  addTiles() %>% 
  addPolygons(data=carreaux_select %>% st_transform(4326),
              label = "Carreaux sélectionnés") %>% 
  addPolygons(data=buffer_white %>% st_transform(4326) ,
              color="red",
              fill=F,
              label = "Zone à façon")

Cette carte représente les carreaux issus de Filosofi carroyé en 2015 qui intersectent notre zone d’intérêt. À noter que ces carreaux sont « penchés » : ceci est la résultant de leur reprojection en Lambert 93. En effet, ces carreaux ont été produits en projection LAEA (standard européen), ce qui leur donne cet aspect une fois reprojetés dans le standard français.

Maintenant que nous avons sélectionné ces carreaux, il suffit de calculer l’indicateur souhaité.

tx_logSoc_white <- sum(carreaux_select$Log_soc)/sum(carreaux_select$Men)
cat("On compte ",
    round(tx_logSoc_white*100),
    "% de logements sociaux dans un rayon de 1km autour du white")

## On compte  33 % de logements sociaux dans un rayon de 1km autour du white

Vous êtes désormais capables de carroyer et lisser toutes les données que vous aurez à votre disposition !

Reproducibilité

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
## [1] leaflet_2.1.0     aws.s3_0.3.21     mapview_2.10.0    mapsf_0.4.0      
## [5] btb_0.1.30.3      sf_1.0-6          dplyr_1.0.8       unilur_0.4.0.9100
## [9] knitr_1.37       
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.2              jsonlite_1.8.0          RcppParallel_5.1.5     
##  [4] assertthat_0.2.1        sp_1.4-6                highr_0.9              
##  [7] stats4_4.1.2            renv_0.15.2             yaml_2.3.5             
## [10] pillar_1.7.0            lattice_0.20-45         glue_1.6.2             
## [13] uuid_1.0-3              digest_0.6.29           RColorBrewer_1.1-2     
## [16] colorspace_2.0-3        leaflet.providers_1.9.0 htmltools_0.5.2        
## [19] pkgconfig_2.0.3         raster_3.5-15           purrr_0.3.4            
## [22] scales_1.1.1            webshot_0.5.2           brew_1.0-7             
## [25] svglite_2.1.0           terra_1.5-21            satellite_1.0.4        
## [28] tibble_3.1.6            proxy_0.4-26            generics_0.1.2         
## [31] farver_2.1.0            ellipsis_0.3.2          cli_3.2.0              
## [34] magrittr_2.0.2          crayon_1.5.0            evaluate_0.15          
## [37] fansi_1.0.2             xml2_1.3.3              class_7.3-19           
## [40] tools_4.1.2             lifecycle_1.0.1         stringr_1.4.0          
## [43] munsell_0.5.0           compiler_4.1.2          jquerylib_0.1.4        
## [46] e1071_1.7-9             systemfonts_1.0.4       rlang_1.0.1            
## [49] classInt_0.4-3          units_0.8-0             grid_4.1.2             
## [52] leafpop_0.1.0           htmlwidgets_1.5.4       aws.signature_0.6.0    
## [55] crosstalk_1.2.0         leafem_0.1.6            base64enc_0.1-3        
## [58] rmarkdown_2.11          codetools_0.2-18        DBI_1.1.2              
## [61] curl_4.3.2              R6_2.5.1                fastmap_1.1.0          
## [64] utf8_1.2.2              KernSmooth_2.23-20      stringi_1.7.6          
## [67] Rcpp_1.0.8              vctrs_0.3.8             png_0.1-7              
## [70] tidyselect_1.1.2        xfun_0.29

Formation analyse urbaine - Travaux pratiques sur R (carroyage, lissage, zones à façon) [TUTORIEL]

Kim Antunez et Julien Pramil

11 mars 2021